## Replication code for
## Weidmann & Ward: Predicting Conflict in Space and Time
## JCR, forthcoming

## This version: March 28, 2010
## Nils B. Weidmann (nils.weidmann@gmail.com)

# Though packaged as such ('streg'), the code to run the spatial-temporal 
# model is not an official R library. Rather, it is available
# in the replication archive and needs to be installed from
# a local source (using R's package manager). It requires the 
# 'rJava' package to be present.
library(streg)
library(ROCR)

# set your working directory here
setwd("<your working dir>")

# load functions
source("functions.R")
source("separationplot.R")

# read data
data <- read.table("bosnia_replication.tab", header=T)
W <-  as.matrix(read.table("bosnia_wmat.txt", header=F, check.names=F))

# make sure the dataset is sorted correctly
data <- streg.sort(data, "month", "MUNID")

# we create the lags
data$CONFLICT.t1 <- streg.tlag(data, "month", "CONFLICT", order=1)
data$CONFLICT.t2 <- streg.tlag(data, "month", "CONFLICT", order=2)
data$CONFLICT.s1 <- streg.slag(data, "month", "CONFLICT", W)
data <- na.omit(data)

## Model estimation

# run stlogit function
model.st <- stlogit(CONFLICT ~ TOTPOP91 + ELF91 + BORDER + terrmean + CONFLICT.t1 + CONFLICT.t2, data=data, unitvar="MUNID", timevar="month", W, seed=2)
model.st

# run temporal logit model
model.t <- glm(CONFLICT ~ TOTPOP91 + ELF91 + BORDER + terrmean + CONFLICT.t1 + CONFLICT.t2, data=data, family=binomial)
summary(model.t)

## In-sample prediction using observed spatial lag

# stlogit
theta <- model.st$theta
eta <- exp(theta[1] + theta[2]*data$TOTPOP91 + theta[3]*data$ELF91 + theta[4]*data$BORDER + theta[5]*data$terrmean + theta[6]*data$CONFLICT.t1 + theta[7]*data$CONFLICT.t2 + theta[8]*data$CONFLICT.s1)
data$phat.st.real <- eta/(1+eta)

# logit
data$phat.t.real <- predict(model.t, type = "response")

# create ROC plot
pred.st.real <- prediction(data$phat.st.real, data$CONFLICT)
perf.st.real <- performance(pred.st.real, measure = "tpr", x.measure = "fpr")
plot(perf.st.real, lty=1, lwd=2)
pred.t.real <- prediction(data$phat.t.real, data$CONFLICT)
perf.t.real <- performance(pred.t.real, measure = "tpr", x.measure = "fpr")
plot(perf.t.real, lty=2, lwd=2, add=T)

# AUC
perf <- performance(pred.st.real, measure = "auc")
perf
perf <- performance(pred.t.real, measure = "auc")
perf

# Separation plots
separationplot(data$phat.t.real, data$CONFLICT, mar=0, lwd=2)
separationplot(data$phat.st.real, data$CONFLICT, mar=0, lwd=2)

## In-sample with simulation
# stlogit
set.seed(1)
for (m in 3:44) {
	cat("Month ", m, "\n")
	sim.month <- simulate.month(data[data$month==m,], model.st$theta, 100, 1000) 
  	data$phat.st.sim[data$month==m] <- sim.month$phats  
  	sse <- sim.month$sse
	#title <- paste("Month", m)
	#pdf(paste(title, ".pdf", sep=""))
	#plot(sse, main=title, xlab="Iteration", ylab="SSE", type="l", lwd=2)
	#dev.off()
}


# ROC
pred <- prediction(data$phat.st.sim, data$CONFLICT)
perf <- performance(pred, measure = "tpr", x.measure = "fpr")
plot(perf, lty=1, lwd=2)
perf <- performance(pred, measure = "auc")
perf
pred <- prediction(data$phat.t.real, data$CONFLICT)
perf <- performance(pred, measure = "tpr", x.measure = "fpr")
plot(perf, lty=2, lwd=2, add=T)
perf <- performance(pred, measure = "auc")
perf 


## Out-of-sample
# logit
for (m in 23:44) {
	train.set <- data[data$month<m,]
	test.set <- data[data$month==m,] 
  	ple <- glm(CONFLICT ~ TOTPOP91 + ELF91 + BORDER + terrmean + CONFLICT.t1 + CONFLICT.t2, data=train.set, family=binomial)
	data$phat.t.os[data$month==m] <- predict(ple, type="response", newdata=test.set) 
}

# spatial temporal
set.seed(1)
for (m in 23:44) {
	cat("Month ", m, "\n")
	train.set <- data[data$month<m,]
	test.set <- data[data$month==m,] 
	m1 <- stlogit(CONFLICT ~ TOTPOP91 + ELF91 + BORDER + terrmean + CONFLICT.t1 + CONFLICT.t2, data=train.set, unitvar="MUNID", timevar="month", W)
	sim.month <- simulate.month(test.set, m1$theta, 100, 1000) 
	data$phat.st.os[data$month==m] <- sim.month$phats
}

# evaluate predictions
pred <- prediction(data$phat.st.os[data$month>=23], data$CONFLICT[data$month>=23])
perf <- performance(pred, measure = "tpr", x.measure = "fpr")
plot(perf, lty=1, lwd=2)
perf <- performance(pred, measure = "auc")
perf
pred <- prediction(data$phat.t.os[data$month>=23], data$CONFLICT[data$month>=23])
perf <- performance(pred, measure = "tpr", x.measure = "fpr")
plot(perf, lty=2, lwd=2, add=T)
perf <- performance(pred, measure = "auc")
perf

# separation plots
separationplot(data$phat.t.os[data$month>=23], data$CONFLICT[data$month>=23], lwd=2)
separationplot(data$phat.st.os[data$month>=23], data$CONFLICT[data$month>=23], lwd=2)




